home *** CD-ROM | disk | FTP | other *** search
/ The CICA Windows Explosion! / The CICA Windows Explosion! - Disc 2.iso / misc / gs261src.zip / gxpcopy.c < prev    next >
C/C++ Source or Header  |  1993-05-13  |  17KB  |  501 lines
  1. /* Copyright (C) 1992, 1993 Aladdin Enterprises.  All rights reserved.
  2.  
  3. This file is part of Ghostscript.
  4.  
  5. Ghostscript is distributed in the hope that it will be useful, but
  6. WITHOUT ANY WARRANTY.  No author or distributor accepts responsibility
  7. to anyone for the consequences of using it or for whether it serves any
  8. particular purpose or works at all, unless he says so in writing.  Refer
  9. to the Ghostscript General Public License for full details.
  10.  
  11. Everyone is granted permission to copy, modify and redistribute
  12. Ghostscript, but only under the conditions described in the Ghostscript
  13. General Public License.  A copy of this license is supposed to have been
  14. given to you along with Ghostscript so you can know your rights and
  15. responsibilities.  It should be in a file named COPYING.  Among other
  16. things, the copyright notice and this notice must be preserved on all
  17. copies.  */
  18.  
  19. /* gxpcopy.c */
  20. /* Path copying and flattening for Ghostscript library */
  21. #include "math_.h"
  22. #include "gx.h"
  23. #include "gserrors.h"
  24. #include "gxfixed.h"
  25. #include "gxarith.h"
  26. #include "gzpath.h"
  27.  
  28. /* Forward declarations */
  29. private int copy_path(P4(const gx_path *, gx_path *, fixed, int));
  30. private int flatten_recur(P8(gx_path *,
  31.   fixed, fixed, fixed, fixed, fixed, fixed, fixed));
  32. private int flatten_sample(P8(gx_path *, int,
  33.   fixed, fixed, fixed, fixed, fixed, fixed));
  34.  
  35. /* Copy a path */
  36. int
  37. gx_path_copy(const gx_path *ppath_old, gx_path *ppath, int init)
  38. {    return copy_path(ppath_old, ppath, max_fixed, init);
  39. }
  40.  
  41. /* Flatten a path. */
  42. /* If flatness is zero, use sampling rather than subdivision: */
  43. /* this is important for Type 1 fonts. */
  44. private fixed
  45. scale_flatness(floatp flatness)
  46. {    fixed scaled_flat = float2fixed(flatness);
  47.     return (scaled_flat > int2fixed(100) ? int2fixed(100) :
  48.         scaled_flat < fixed_half ? fixed_0 :
  49.         scaled_flat);
  50. }
  51. int
  52. gx_path_flatten(const gx_path *ppath_old, gx_path *ppath, floatp flatness, int in_BuildCharGlyph)
  53. {    return copy_path(ppath_old, ppath,
  54.              (in_BuildCharGlyph ? fixed_0 : scale_flatness(flatness)),
  55.              1);
  56. }
  57.  
  58. /* Add a flattened curve to a path. */
  59. int
  60. gx_path_add_flattened_curve(gx_path *ppath,
  61.   fixed x1, fixed y1, fixed x2, fixed y2, fixed x3, fixed y3,
  62.   floatp flatness)
  63. {    return flatten_recur(ppath, x1, y1, x2, y2, x3, y3,
  64.                  scale_flatness(flatness));
  65. }
  66.  
  67. /* Copy a path, optionally flattening it. */
  68. /* If the copy fails, free the new path. */
  69. private int
  70. copy_path(const gx_path *ppath_old, gx_path *ppath, fixed scaled_flat,
  71.   int init)
  72. {    gx_path old;
  73.     const segment *pseg;
  74.     int code;
  75. #ifdef DEBUG
  76. if ( gs_debug['p'] )
  77.     gx_dump_path(ppath_old, "before copy_path");
  78. #endif
  79.     old = *ppath_old;
  80.     if ( init )
  81.         gx_path_init(ppath, ppath_old->memory_procs);
  82.     pseg = (const segment *)(old.first_subpath);
  83.     while ( pseg )
  84.        {    switch ( pseg->type )
  85.            {
  86.         case s_start:
  87.             code = gx_path_add_point(ppath, pseg->pt.x, pseg->pt.y);
  88.             break;
  89.         case s_curve:
  90.            {    curve_segment *pc = (curve_segment *)pseg;
  91.             if ( scaled_flat == max_fixed )    /* don't flatten */
  92.                 code = gx_path_add_curve(ppath,
  93.                     pc->p1.x, pc->p1.y,
  94.                     pc->p2.x, pc->p2.y,
  95.                     pc->pt.x, pc->pt.y);
  96.             else
  97.                 code = flatten_recur(ppath,
  98.                     pc->p1.x, pc->p1.y,
  99.                     pc->p2.x, pc->p2.y,
  100.                     pc->pt.x, pc->pt.y,
  101.                     scaled_flat);
  102.             break;
  103.            }
  104.         case s_line:
  105.             code = gx_path_add_line(ppath, pseg->pt.x, pseg->pt.y);
  106.             break;
  107.         case s_line_close:
  108.             code = gx_path_close_subpath(ppath);
  109.             break;
  110.            }
  111.         if ( code )
  112.            {    gx_path_release(ppath);
  113.             if ( ppath == ppath_old )
  114.                 *ppath = old;
  115.             return code;
  116.            }
  117.         pseg = pseg->next;
  118.     }
  119.     if ( !old.subpath_open && old.position_valid )
  120.         gx_path_add_point(ppath, old.position.x, old.position.y);
  121. #ifdef DEBUG
  122. if ( gs_debug['p'] )
  123.     gx_dump_path(ppath, "after copy_path");
  124. #endif
  125.     return 0;
  126. }
  127. /* Internal routine to flatten a curve. */
  128. /* This calls itself recursively, using binary subdivision, */
  129. /* until the approximation is good enough to satisfy the */
  130. /* flatness requirement.  The starting point is ppath->position, */
  131. /* which gets updated as line segments are added. */
  132.  
  133. /* Maximum number of points for sampling if we want accurate rasterizing. */
  134. /* (num_sample_max - 1)^3 must fit into an int with a bit to spare. */
  135. #if arch_sizeof_int == 2
  136. #  define num_sample_max 26
  137. #else
  138. #  define num_sample_max (1 << ((size_of(int) * 8 - 2) / 3))
  139. #endif
  140.  
  141. /* Table of f(i) = 256 * sqrt(1 + (i/64)^2). */
  142. /* This is good to within 1% or better. */
  143. #define sqrt_index_shift 6        /* scaling of index */
  144. #define sqrt_value_shift 8        /* scaling of value */
  145. private int scaled_sqrt_tab[65] =
  146.    {    256, 256, 256, 256, 256, 256, 257, 257,
  147.     257, 258, 259, 259, 260, 261, 262, 262,
  148.     263, 264, 265, 267, 268, 269, 270, 272,
  149.     273, 274, 276, 277, 279, 281, 282, 284,
  150.     286, 288, 289, 291, 293, 295, 297, 299,
  151.     301, 304, 306, 308, 310, 312, 315, 317,
  152.     320, 322, 324, 327, 329, 332, 334, 337,
  153.     340, 342, 345, 348, 350, 353, 356, 359,
  154.     362
  155.    };
  156.  
  157. private int
  158. flatten_recur(gx_path *ppath,
  159.   fixed x1, fixed y1, fixed x2, fixed y2, fixed x3, fixed y3,
  160.   fixed scaled_flat)
  161. {    fixed
  162.       x0 = ppath->position.x,
  163.       y0 = ppath->position.y;
  164. top:
  165. #ifdef DEBUG
  166. if ( gs_debug['2'] )
  167.     dprintf4("[2]x0=%f y0=%f x1=%f y1=%f\n",
  168.          fixed2float(x0), fixed2float(y0),
  169.          fixed2float(x1), fixed2float(y1)),
  170.     dprintf4("   x2=%f y2=%f x3=%f y3=%f\n",
  171.          fixed2float(x2), fixed2float(y2),
  172.          fixed2float(x3), fixed2float(y3));
  173. #endif
  174.     /*
  175.      * Compute the maximum distance of the curve from
  176.      * the line (x0,y0)->(x3,y3).  We do this conservatively
  177.      * by observing that the curve is enclosed by the
  178.      * quadrilateral of its control points, so we simply
  179.      * compute the distances of (x1,y1) and (x2,y2)
  180.      * from the line.  Letting dx = x3-x0 and dy = y3-y0,
  181.      * the distance of (xp,yp) from the line is
  182.      * abs(N)/sqrt(D), where N = dy*(xp-x0)-dx*(yp-y0) and
  183.      * D = dx*dx+dy*dy; hence we want to test abs(N) <= sqrt(D)*F,
  184.      * where F is the flatness parameter from the graphics state.
  185.      * We can do this more efficiently by letting t=dy/dx, and
  186.      * testing abs(N1) <= sqrt(D1)*f, where N1=t*(xp-x0)-(yp-y0) and
  187.      * D1 = 1+t*t.  If dx < dy, we swap x and y for this
  188.      * computation.  This guarantees abs(t) <= 1, which allows us to
  189.      * compute sqrt(1+t*t) by table lookup on the high bits of abs(t).
  190.      *
  191.      * To avoid replacing shallow curves by short straight lines,
  192.      * we halve the flatness if the curve is very short.
  193.      * We don't do anything about long, nearly flat curves.
  194.      *
  195.      * Note that if scaled_flat is very small, we don't do any of this;
  196.      * instead, we just check whether abs(dx) and abs(dy) are
  197.      * small enough to switch over to sampling rather than dividing.
  198.       */
  199.      { fixed dx3 = x3 - x0;
  200.        fixed adx3 = any_abs(dx3);
  201.        fixed dy3 = y3 - y0;
  202.        fixed ady3 = any_abs(dy3);
  203.        fixed flat = scaled_flat;
  204.        /* We have to be quite careful to ensure that */
  205.        /* none of the multiplications will overflow. */
  206. #define short_max 0x7ff0L
  207. #define reduce_3(ad3, maxv)\
  208.   while ( ad3 > maxv )\
  209.     adx3 >>= 1, ady3 >>= 1,\
  210.     dx3 = arith_rshift_1(dx3), dy3 = arith_rshift_1(dy3)
  211. #define reduce_d(d)\
  212.   for ( shift = 0; (d < 0 ? d < -short_max : d > short_max); shift++ )\
  213.     d = arith_rshift_1(d)
  214.        if ( adx3 > ady3 )
  215.         {    fixed d, dx, dy, dist;
  216.         int shift;
  217.         if ( scaled_flat == 0 )
  218.         { if ( adx3 < int2fixed(num_sample_max - 1) / 2 )
  219.           { int n = fixed2int_var_rounded(adx3) * 2 + 1;
  220.             int code = flatten_sample(ppath, max(n, 3), x1, y1,
  221.                           x2, y2, x3, y3);
  222.             if ( code <= 0 ) return code;
  223.           }
  224.           goto sub;
  225.         }
  226.         if ( adx3 < 16 ) flat >>= 1;
  227.         reduce_3(ady3, short_max);
  228.         d = (scaled_sqrt_tab[(ady3 << sqrt_index_shift) / adx3] * flat) >> sqrt_value_shift;
  229.         dx = x1 - x0, dy = y1 - y0;
  230.         reduce_d(dx);
  231.         if ( ((dist = ((dx * dy3 / dx3) << shift) - dy) < 0 ?
  232.               -dist : dist) > d )
  233.           goto sub;    /* not flat enough */
  234.         dx = x2 - x0, dy = y2 - y0;
  235.         reduce_d(dx);
  236.         if ( ((dist = ((dx * dy3 / dx3) << shift) - dy) < 0 ?
  237.               -dist : dist) > d )
  238.           goto sub;    /* not flat enough */
  239.         }
  240.        else if ( ady3 != 0 )
  241.         {    fixed d, dy, dx, dist;
  242.         int shift;
  243.         if ( scaled_flat == 0 )
  244.         { if ( ady3 < int2fixed(num_sample_max - 1) / 2 )
  245.           { int n = fixed2int_var_rounded(ady3) * 2 + 1;
  246.             int code = flatten_sample(ppath, max(n, 3), x1, y1,
  247.                           x2, y2, x3, y3);
  248.             if ( code <= 0 ) return code;
  249.           }
  250.           else goto sub;
  251.         }
  252.         if ( ady3 < 16 ) flat >>= 1;
  253.         reduce_3(adx3, short_max);
  254.         d = (scaled_sqrt_tab[(adx3 << sqrt_index_shift) / ady3] * flat) >> sqrt_value_shift;
  255.         dy = y1 - y0, dx = x1 - x0;
  256.         reduce_d(dy);
  257.         if ( ((dist = ((dy * dx3 / dy3) << shift) - dx) < 0 ?
  258.               -dist : dist) > d )
  259.           goto sub;    /* not flat enough */
  260.         dy = y2 - y0, dx = x2 - x0;
  261.         reduce_d(dy);
  262.         if ( ((dist = ((dy * dx3 / dy3) << shift) - dx) < 0 ?
  263.               -dist : dist) > d )
  264.           goto sub;    /* not flat enough */
  265.         }
  266.        else                /* adx3 = ady3 = 0 */
  267.         {    /* (x0,y0) is the same point as (x3,y3). */
  268.         /* This is an anomalous case.  If the entire curve */
  269.         /* is a single point, stop now, otherwise subdivide. */
  270.         if ( x1 != x0 || y1 != y0 || x2 != x0 || y2 != y0 )
  271.           goto sub;
  272.         }
  273.      }
  274.     /* Curve is flat enough.  Add a line and exit. */
  275.     if_debug2('2', "[2]\t*** x=%f, y=%f ***\n",
  276.           fixed2float(x3), fixed2float(y3));
  277.     return gx_path_add_line(ppath, x3, y3);
  278.  
  279.     /* Curve isn't flat enough.  Break into two pieces and recur. */
  280.     /* Algorithm is from "The Beta2-split: A special case of the */
  281.     /* Beta-spline Curve and Surface Representation," B. A. Barsky */
  282.     /* and A. D. DeRose, IEEE, 1985, courtesy of Crispin Goswell. */
  283. sub:
  284.     /* We have to define midpoint carefully to avoid overflow. */
  285.     /* (If it overflows, something really pathological is going on, */
  286.     /* but we could get infinite recursion that way....) */
  287. #define midpoint(a,b)\
  288.   (arith_rshift_1(a) + arith_rshift_1(b) + ((a) & (b) & 1))
  289.    {    fixed x01 = midpoint(x0, x1), y01 = midpoint(y0, y1);
  290.     fixed x12 = midpoint(x1, x2), y12 = midpoint(y1, y2);
  291.     fixed x02 = midpoint(x01, x12), y02 = midpoint(y01, y12);
  292.     int code;
  293.     /* Update x/y1, x/y2, and x/y0 now for the second half. */
  294.     x2 = midpoint(x2, x3), y2 = midpoint(y2, y3);
  295.     x1 = midpoint(x12, x2), y1 = midpoint(y12, y2);
  296.     code = flatten_recur(ppath, x01, y01, x02, y02,
  297.         (x0 = midpoint(x02, x1)), (y0 = midpoint(y02, y1)),
  298.         scaled_flat);
  299.     if ( code < 0 ) return code;
  300.    }    goto top;
  301. }
  302.  
  303. /* Flatten a segment of the path by repeated sampling. */
  304. /* n is the number of points to sample, including the endpoints; */
  305. /* we require n >= 3. */
  306. /* If n or any of the coefficient values are too large, */
  307. /* return 1 and let the recursive subdivider whittle them down further. */
  308. private int
  309. flatten_sample(gx_path *ppath, int n,
  310.   fixed x1, fixed y1, fixed x2, fixed y2, fixed x3, fixed y3)
  311. {    const fixed
  312.         x0 = ppath->position.x,
  313.         y0 = ppath->position.y;
  314.     /* We spell out some multiplies by 3, */
  315.     /* for the benefit of compilers that don't optimize this. */
  316.     const fixed
  317.         x01 = x1 - x0,
  318.         cx = (x01 << 1) + x01,        /* 3*(x1-x0) */
  319.         x12 = x2 - x1,
  320.         bx = (x12 << 1) + x12 - cx,    /* 3*(x2-2*x1+x0) */
  321.         ax = x3 - bx - cx - x0;        /* x3-3*x2+3*x1-x0 */
  322.     const fixed
  323.         y01 = y1 - y0,
  324.         cy = (y01 << 1) + y01,
  325.         y12 = y2 - y1,
  326.         by = (y12 << 1) + y12 - cy,
  327.         ay = y3 - by - cy - y0;
  328.     const int
  329.         n1 = n - 1,
  330.         n12 = n1 * n1,
  331.         n13 = n12 * n1;
  332.     fixed ptx = x0, pty = y0;
  333.     fixed x, y;
  334.     /*
  335.      * If all the coefficients lie between min_fast and max_fast,
  336.      * we can do everything in fixed point.  In this case we compute
  337.      * successive values by finite differences, using the formulas:
  338.         x(t) =
  339.           a*t^3 + b*t^2 + c*t + d =>
  340.         dx(t) = x(t+e)-x(t) =
  341.           a*(3*t^2*e + 3*t*e^2 + e^3) + b*(2*t*e + e^2) + c*e =
  342.           (3*a*e)*t^2 + (3*a*e^2 + 2*b*e)*t + (a*e^3 + b*e^2 + c*e) =>
  343.         d2x(t) = dx(t+e)-dx(t) =
  344.           (3*a*e)*(2*t*e + e^2) + (3*a*e^2 + 2*b*e)*e =
  345.           (6*a*e^2)*t + (6*a*e^3 + 2*b*e^2) =>
  346.         d3x(t) = d2x(t+e)-d2x(t) =
  347.           6*a*e^3;
  348.         x(0) = d, dx(0) = (a*e^3 + b*e^2 + c*e),
  349.           d2x(0) = 6*a*e^3 + 2*b*e^2;
  350.      * In these formulas, e = 1/n1; of course, there are separate
  351.      * computations for the x and y values.
  352.      */
  353. #define max_fast (max_fixed / 6)
  354. #define min_fast (-max_fast)
  355.     int i;
  356.     /*
  357.      * We do exact rational arithmetic to avoid accumulating error.
  358.      * Each quantity is represented as I+R/n13, where I is an "integer"
  359.      * and the "remainder" R lies in the range 0 <= R < n13.  Note that
  360.      * R may temporarily exceed n13; for this reason, we require that
  361.      * n13 have at least one free high-order bit (below the sign).
  362.      */
  363.     fixed idx, idy, id2x, id2y, id3x, id3y;        /* I */
  364.     int rx, ry, rdx, rdy, rd2x, rd2y, rd3x, rd3y;
  365.  
  366.     if_debug6('2', "[2]ax=%f bx=%f cx=%f\n   ay=%f by=%f cy=%f\n",
  367.           fixed2float(ax), fixed2float(bx), fixed2float(cx),
  368.           fixed2float(ay), fixed2float(by), fixed2float(cy));
  369. #define in_range(v) (v < max_fast && v > min_fast)
  370.     if ( n1 < num_sample_max &&    /* so n13 fits */
  371.          in_range(ax) && in_range(bx) && in_range(cx) &&
  372.          in_range(ay) && in_range(by) && in_range(cy)
  373.        )
  374.     {    x = x0, y = y0;
  375.         rx = ry = 0;
  376.         /* Fast check for n == 3, a common special case */
  377.         /* for small characters. */
  378.         if ( n == 3 )
  379. #define poly2(a,b,c)\
  380.   arith_rshift_1(arith_rshift_1(arith_rshift_1(a) + b) + c)
  381.             idx = poly2(ax, bx, cx),
  382.             idy = poly2(ay, by, cy),
  383.             rdx = rdy = 0;
  384. #undef poly2
  385.         else
  386.         {    fixed bx2 = bx << 1, by2 = by << 1;
  387.             fixed ax6 = ((ax << 1) + ax) << 1,
  388.                   ay6 = ((ay << 1) + ay) << 1;
  389.             fixed qx, qy;
  390. #define adjust_rem(r, q)\
  391.   if ( r < 0 ) r += n13, q--
  392. #define adjust_rem_loop(r, q)\
  393.   while ( r < 0 ) r += n13, q--;\
  394.   while ( r >= n13 ) r -= n13, q++
  395.             /* We can compute all the remainders as ints, */
  396.             /* because we know they are less than n13. */
  397.             /* bx/y terms */
  398.             id2x = bx2 / n12, id2y = by2 / n12;
  399.             rd2x = ((int)bx2 - (int)id2x * n12) * n1,
  400.               rd2y = ((int)by2 - (int)id2y * n12) * n1;
  401.             idx = bx / n12, idy = by / n12;
  402.             rdx = ((int)bx - (int)idx * n12) * n1,
  403.               rdy = ((int)by - (int)idy * n12) * n1;
  404.             /* cx/y terms */
  405.             idx += qx = cx / n1, idy += qy = cy / n1;
  406.             rdx += ((int)cx - (int)qx * n1) * n12,
  407.               rdy += ((int)cy - (int)qy * n1) * n12;
  408.             /* ax/y terms */
  409.             idx += qx = ax / n13, idy += qy = ay / n13;
  410.             rdx += (int)ax - (int)qx * n13,
  411.               rdy += (int)ay - (int)qy * n13;
  412.             id2x += id3x = ax6 / n13, id2y += id3y = ay6 / n13;
  413.             rd2x += rd3x = (int)ax6 - (int)id3x * n13,
  414.               rd2y += rd3y = (int)ay6 - (int)id3y * n13;
  415.             adjust_rem_loop(rdx, idx);
  416.             adjust_rem_loop(rdy, idy);
  417.             adjust_rem_loop(rd2x, id2x);
  418.             adjust_rem_loop(rd2y, id2y);
  419.             adjust_rem(rd3x, id3x);
  420.             adjust_rem(rd3y, id3y);
  421. #undef adjust_rem
  422. #undef adjust_rem_loop
  423.         }
  424.     }
  425.     else
  426.     {    /* Subdivide some more. */
  427.         return 1;
  428.     }
  429.     if_debug1('2', "[2]sampling n=%d\n", n);
  430.     for ( i = 1; ; )
  431.     {    int code;
  432. #ifdef DEBUG
  433. if ( gs_debug['2'] )
  434.         dprintf4("[2]dx=%f+%d, dy=%f+%d\n",
  435.              fixed2float(idx), rdx,
  436.              fixed2float(idy), rdy),
  437.         dprintf4("   d2x=%f+%d, d2y=%f+%d\n",
  438.              fixed2float(id2x), rd2x,
  439.              fixed2float(id2y), rd2y),
  440.         dprintf4("   d3x=%f+%d, d3y=%f+%d\n",
  441.              fixed2float(id3x), rd3x,
  442.              fixed2float(id3y), rd3y);
  443. #endif
  444. #define accum(i, r, di, dr)\
  445.   if ( (r += dr) >= n13 ) r -= n13, i += di + 1;\
  446.   else i += di
  447.         accum(x, rx, idx, rdx);
  448.         accum(y, ry, idy, rdy);
  449.         if_debug3('2', "[2]%s x=%g, y=%g\n",
  450.               (((x ^ ptx) | (y ^ pty)) & float2fixed(-0.5) ?
  451.                "add" : "skip"),
  452.               fixed2float(x), fixed2float(y));
  453.         /* Skip very short segments */
  454.         if ( ((x ^ ptx) | (y ^ pty)) & float2fixed(-0.5) )
  455.         {    if ( (code = gx_path_add_line(ppath, x, y)) < 0 )
  456.                 return code;
  457.             ptx = x, pty = y;
  458.         }
  459.         if ( ++i >= n1 )
  460.             break;        /* don't bother with last accum */
  461.         accum(idx, rdx, id2x, rd2x);
  462.         accum(id2x, rd2x, id3x, rd3x);
  463.         accum(idy, rdy, id2y, rd2y);
  464.         accum(id2y, rd2y, id3y, rd3y);
  465. #undef accum
  466.     }
  467.     if_debug2('2', "[2]last x=%g, y=%g\n",
  468.           fixed2float(x3), fixed2float(y3));
  469.     return gx_path_add_line(ppath, x3, y3);
  470. }
  471.  
  472. /*
  473.  *    The rest of this file is an analysis that will eventually
  474.  *    allow us to rasterize curves on the fly, by finding points
  475.  *    where Y reaches a local maximum or minimum, which allows us to
  476.  *    divide the curve into locally Y-monotonic sections.
  477.  */
  478.  
  479. /*
  480.     Let y(t) = a*t^3 + b*t^2 + c*t + d, 0 <= t <= 1.
  481.     Then dy(t) = 3*a*t^2 + 2*b*t + c.
  482.     y(t) has a local minimum or maximum (or inflection point)
  483.     precisely where dy(t) = 0.  Now the roots of dy(t) are
  484.         ( -2*b +/- sqrt(4*b^2 - 12*a*c) ) / 6*a
  485.        =    ( -b +/- sqrt(b*2 - 3*a*c) ) / 3*a
  486.     (Note that real roots exist iff b^2 >= 3*a*c.)
  487.     We want to know if these lie in the range (0..1).
  488.     (We don't care about the endpoints.)  Call such a root
  489.     a "valid zero."  We proceed as follows:
  490.         If sign(3*a + 2*b + c) ~= sign(c), a valid zero exists.
  491.         If sign(a) = sign(b), no valid zero exists.
  492.     Otherwise, we look for a local extremum of dy(t) by observing
  493.         d2y(t) = 6*a*t + 2*b
  494.     which has a zero only at
  495.         t1 = -b / 3*a
  496.     Now if t1 <= 0 or t1 >= 1, no valid zero exists.  Otherwise,
  497.     we compute
  498.         dy(t1) = c - (b^2 / 3*a)
  499.     Then a valid zero exists (at t1) iff sign(dy(t1)) ~= sign(c).
  500.  */
  501.